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Abstract 

This paper discusses the mathematical framework for designing methods of large defor- 
mation matching (LDM) for image registration in computational anatomy. After reviewing 
the geometrical framework of LDM image registration methods, a theorem is proved show- 
ing that these methods may be designed by using the actions of diffeomorphisms on the 
image data structure to define their associated momentum representations as (cotangent 
lift) momentum maps. To illustrate its use, the momentum map theorem is shown to re- 
cover the known algorithms for matching landmarks, scalar images and vector fields. After 
briefly discussing the use of this approach for Diffusion Tensor (DT) images, we explain 
how to use momentum maps in the design of registration algorithms for more general data 
structures. For example, we extend our methods to determine the corresponding momen- 
tum map for registration using semidirect product groups, for the purpose of matching 
images at two different length scales. Finally, we discuss the use of momentum maps in 
the design of image registration algorithms when the image data is defined on manifolds 
instead of vector spaces. 
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1 Introduction 

Large deformation diffeomorphic matching methods (LDM) for image registration are based on 
minimizing the sum of a kinetic energy metric, plus a penalty term. The former ensures that 
the deformation follows an optimal path, while the latter ensures an acceptable tolerance in 
image mismatch. The LDM approaches were introduced and systematically developed in Trouve 
[29, 30], Dupuis et al. [10], Joshi and Miller [19], Miller and Younes [22], Beg [5], and Beg et al. [7]. 
See Miller et al. [23] for an extensive review of this development. The LDM approach fits within 
Grenander's [12] deformable template paradigm for image registration. Grenander's paradigm, 
in turn, is a development of the biometric strategy introduced by Thompson [28] of comparing 
a template image Jq to a target image Ii by finding a smooth transformation that maps the 
template to the target. This transformation is assumed to belong to a Lie group G that acts on 
the set of images V containing Iq and Ii. The effect of the transformation on the data structure 
is called the action G x V V of the Lie group G on the set V. For example, the action of 
g E G on Iq E V is denoted as glo G V. 

The objective of LDM is not just to determine a deformation gi E G such that the group 
action gil^ oi gi E G on the template Iq E V approximates the target Ii E V to within a 
certain tolerance. Rather, the objective of LDM is to find the optimal path gt E G continuously 
parametrized by time t G M that smoothly deforms Jq through If = gtl^ to gil^. The optimal 
path gt E G is defined as the path that costs the least in time-integrated kinetic energy for a 
given tolerance. Hence, the deformable template method may be formulated as an optimization 
problem based on a trade-off between the following two properties: (i) the tolerance for inexact 
matching between the final deformed template g^lQ and the target template Ii, and (2) the cost of 
time-integrated kinetic energy of the rate of deformation along the path gt- The former is defined 
by assigning a norm || ■ || : V — R to measure the mismatch \\gtIo — /i|| between the two images. 
The latter is obtained by choosing a Riemannian metric | ■ | : TG M that defines the kinetic 
energy on the tangent space TG of the group G. In this setting, a notion of distance between 
two images emerges, that allows one to compare similarity of images in terms of transformations. 
This is the setting for the development of computational anatomy using the inexact template 
matching approach for the registration of images. For more details and background about LDM, 
see Miller and Younes [22], Miller et al. [23], Beg [5], and Beg et al. [7]. 

In applications of LDM to the analysis of features in bio-medical images, the optimal path 
gt is naturally chosen from among the diffeomorphic transformations G = Diff(r2) of an open, 
bounded domain Q. The domain Q will be taken to be the ambient space in which the anatomy 
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is located. Recall that a diffeomorphism g G Diff(r2) is a smooth invertible map (i.e., a invertible 
function that maps the domain Q onto itself) whose inverse is also smooth. The one-to-one 
property of these transformations ensures that disjoint sets remain disjoint, so that, e.g., no 
fusion of points occurs under LDM. Continuity of the diffeomorphisms ensures that connected sets 
remain connected. Smoothness of these transformations ensures preservation of the smoothness of 
boundaries of the anatomical objects in bio-medical images. The invertibility of diffeomorphisms 
and their stability under composition also allows one to regard Diff(r2) formally as a Lie group. 

Different types of bio-medical images contain various types of information that may be rep- 
resented in a number of geometrically different types of data structures. For example, the data 
structures for MR images are scalar functions, or densities, while data obtained DT-MRI can be 
represented as symmetric tensor fields. Naturally, the design of image registration algorithms 
based on the theory of transformations must take differences in data structure into account. 

Registration of DT-MRI data - necessary for the quantitative analysis of anatomical features 
such as tissue geometry and local fiber orientation - is much more complicated than registration 
of scalar image data. This complication arises because local fiber orientation changes under 
a diffeomorphic transformation and this reorientation has to be included properly in the de- 
sign of LDM matching algorithms for DT-MRI. A further complication arises because it is not 
entirely understood how macroscopic deformation influences microscopic properties such as fiber- 
orientation and diffusivity of water. Though significant efforts have been directed at scalar image 
registration, little work has been done on matching tensor images using LDM. For the pioneering 
efforts in the use of LDM with DT-MRI see Alexander et al. [2, 3], Cao et al. [8, 9]. 

In summary, the LDM approach models computational anatomy as a deformation of an initial 
template configuration. The images describing the anatomy are defined on an open bounded set 
Q and the path from the template image Jq to the target image Ji is viewed as the continuous 
deformation It := gtio under the path of diffeomorphic transformations gt G Diff (fi) acting on the 
initial template Iq. Importantly, the optimal path of diffeomorphic transformations gt depends 
on three main factors: namely, how the action gtIo is defined, as well as the definitions of the 
kinetic energy and the tolerance norm. Images representing different types of information may 
transform differently under G = DiS{Q). Hence, the optimal path gt G Diff(M) sought in the 
LDM approach will depend on the geometrical properties of the data structures that represent 
the information in the various types of images. 

In the geometrical framework for the LDM approach, the optimal transformation path gt G 
Diff(M) may be estimated by using the variational optimization method developed in Beg [5] 
and Beg et al. [7]. Namely, the optimal path for the matching diffeomorphism in this problem 
may be obtained from a gradient-descent algorithm based on the directional derivative of the 
cost functional. The cost functional must balance the energy of the deformation path versus 
the tolerance of mismatch, while taking proper account of the transformation properties of the 
image data structure. Other promising methods besides LDM exist, such as the metamorphosis 
approach discussed in Miller and Younes [22], Trouve and Younes [31] and Holm et al. [18]. 
Metamorphosis is a variant of LDM, that allows the evolution It of the image template to deviate 
from pure deformation. It is also a promising method in the LDM family, but its discussion is 
beyond our present scope. 

Our aim in this paper is to show that a simple and universal property of transformation theory, 
called the momentum map can be used to identify and derive the LDM algorithm corresponding 
to any data structure on which diffeomorphisms may act. That is, the momentum map approach 



The momentum map representation of images 



December 14, 2009 4 



enables one to tailor the LDM algorithm to the transformation properties of the data structure 
of the images to be matched. For basic introductions to the momentum map in geometric 
mechanics, see Holm [14] or Marsden and Ratiu [20] . For more extensive treatments see Abraham 
and Marsden [1], Ortega and Ratiu [25]. 

Our interests here focus mainly on deriving the momentum maps corresponding to the various 
types of data structures, rather than developing the matching dynamics that they subsequently 
produce. In particular, we shall discuss how one uses the momentum map approach to cope with 
different data structures, such as densities, vector fields or tensor fields, by recognizing their 
shared properties in a unified geometrical framework. 

Plan of the paper. In Section 2 we begin by discussing the geometry underlying the standard 
algorithm for LDM introduced in Beg et al. [7]. With this motivation we then introduce an 
abstract framework in which to model registration problems. We derive the equivalent of Beg's 
formula in the abstract framework in Theorem 2.5 and show that it has the structure of a 
momentum map. The end of the section is devoted to a discussion of the EPDiff equation and 
the importance of the initial momentum. 

After presenting the abstract framework we apply it in Section 3 to a range of examples 
commonly encountered in computational anatomy: landmarks, scalar images, vector fields and 
symmetric tensor fields arising from DT-MRIs. We emphasize the momentum maps in these 
examples as the main ingredient in our framework and show how to recover results found in the 
literature. 

Section 4 is devoted to a generalization of standard LDM in a different direction. This 
section takes into account the presence of two different length scales in the image and formulates 
a version of LDM that uses a semidirect product of two diffeomorphism groups — one for each 
length scale — to perform the registration. We show that for images defined by scalar functions 
this approach yields a momentum map that is very similar to Beg's formula, except that we use 
the sum of two kernels, instead of only one kernel. 

Besides the formulation of LDM as in Beg et al. [7], other penalty terms have been proposed 
by Beg and Khan [6], Avants et al. [4] and Hart et al. [13]. We show in Section 5 that these 
other proposed penalty terms result in momentum map structures that are similar to those in 
the formulation of Theorem 2.5. 

Our approach can also be generalized to include data structures defined on manifolds that 
do not possess a linear structure. In Section 6 we consider the extensions of the theory required 
to deal with data structure defined on manifolds and apply these extensions in examples. 

2 Geometry of Registration 
2.1 Motivation 

The optimal solution to a non-rigid template matching problem is defined as the shortest, or 
least expensive, path of continuous deformations of one geometric object (template) into another 
one (target). The goal is the find the path of deformations of the template that is shortest, or 
costs the least, for a given tolerance in matching the target. The approach focuses its attention 
on the properties of the action of a Lie group G of transformations on the set of deformable 
templates. The attribution of a cost to this process is based on metrics defined on the tangent 
space TG of the group G, following Grenander's [12] principles. 
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Formulation of LDM 

In the LDM framework this template matching procedure for image registration is formulated 
as follows. Suppose an image, say a medical image, is acquired using MRI, CT, or some other 
imaging technique. To begin, consider the case that the information in an image can be repre- 
sented as a function / : ^2 ^ M, where Q C M*^ is the domain of the image. We denote the data 
structure by writing I E V = J^{^1), the space of smooth functions encoding the information in 
the images. One usually deals with planar [d = 2) or volumetric [d = 3) images. Consider the 
comparison of two images, consisting of a function Jq representing the template image and Ji the 
target image. The goal is to find a transformation : f2 — f2, such that the transformed image 
/q o matches the target image Ji with minimal error, as measured by, say, the norm of 
their difference 

E2(/o,/l) = ||/oO 0-1-/1 Hi.. 

For this purpose, one introduces a time-indexed deformation process, that starts at time t = 
with the template (denoted Iq), and reaches the target Ii at time t = 1. At a given time t during 
this process, the current object It is assumed to be the image of the template, Jq, obtained 
through a sequence of deformations. 

We also want the time-indexed transformation to be regular. To ensure its regularity, we 
require the transformation to be generated as the flow of a smooth time dependent vector fleld 
u : [0,1] X Q ^ Q, i.e. = 0i with 

dt(t)t = ut o (t)t, 0o(x)=x. (2.1) 

We measure the regularity of ut via a kinetic-energy like term 

Ei{ut) = / lutlffdi 
Jo 

where lutln is a norm on the space of vector flelds on Q deflned in terms of a positive self-adjoint 
differential operator L by 

\ut\l^ = {u,Lu)l^ . 

The operator L is commonly chosen as Lu = u — a'^Au. We denote by H this space of vector 
fields. 

Following Beg et al. [7] we can cast the problem of registering /q to Ii as a variational problem. 
Namely, we seek to minimize the cost 

1 

E{ut)= / \ut\^^dt +— \\Ioo^-' -hWl, (2.2) 
Jo 

over all time-dependent vector fields Ut- The transformation 0i is related to the vector field Ut 
via (2.1). A necessary condition for a vector field Ut to be minimal is that the derivative of the 
cost functional E vanishes at Ut, that is DE{ut) = 0. It is shown in Beg et al. [7, theorem 2.1.] 
and Miller et al. [23, theorem 4.1] that DE{ut) = is equivalent to 

Lut = ^\detD<p-l\{J^ - Jl)VJ^, (2.3) 

where 0f^s = (f>t° 07^ and J° = /q ° 0^0' -^1 = h ° 0^1 • This condition is then used in Beg et al. 
[7] to devise a gradient descent algorithm for numerically computing the optimal transformation 

01. 
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Geometric reformulation of LDM 

Formula (2.3) can be reformulated equivalently in a way that emphasizes its geometric nature. 
As we will show in Section 2.2, formula (2.3) is equivalent to 

Lut = {(Pt ■ /o) o (</.i,i ■ {<P, ■ lo - hf) . (2.4) 

This formula can be understood as follows: the first factor (pt'Io is the action of the transformation 
4>t on the image Iq E V = J-'{Q). This is defined as the composition of functions, (pt-h = -^o°0<~^- 
The fiat-operator \) : V V* maps images in V to the objects in V* dual to scalar functions, 
using the inner product on V. (These dual objects are the scalar densities.) To describe such an 
operator, one first needs to choose a convenient space V* in nondegenerate duality with V. We 
choose to identify V* with functions in by using the L^-pairing 

(/,/) := / f{x)Iix)dx, 
Jn 

where dx is a fixed volume element on Q. With this choice, the fiat operator ( b ) is simply the 
identity map on functions. However, it is important that we conceptually distinguish between 
elements in V and in its dual V*. Indeed, the action of a transformation (p on an element in V* 
is the dual action, and does not coincide with the action on V in general. In our example, the 
action on f E V* is 

^■f = \detDr'\ifor')- (2.5) 
To see how this action arises, we need the abstract definition of a dual action, which is 

(0-/,/) = (/,r'-/)- 

The inverse is necessary to ensure that we have a left action: 

0-(^-/) = (0oV;)-/. 
Using this definition and the change of variables formula we see that 

(0-/,/) = (/,r'-/)= [{Io4>)fdx= [ I{fo4>-')\detD(p-'\dx 

Jn Jn 

= (|detD0-i|(/o0-i),/) . 

Therefore, in the second factor 0f i ■ (0i ■ Iq — l{f of equation (2.4), the term (0i ■ Jq — l\f is 
interpreted as a function in V* . Consequently, the action is the dual action given by 

■ (01 ■ Jo - /i)' = |detD0-,i|(J° - J]) . 

It remains to explain the last ingredient, the diamond map 

o:VxV*-^n*. (2.6) 

This is the cotangent-lift momentum map associated to the given action of G on V. Such 
momentum maps are familiar in geometric mechanics; see, e.g.. Holm [14] or Marsden and Ratiu 
[20]. The momentum map (2.6) takes elements ofVx V*, regarded as the cotangent bundle T*V 
of the space of images V, to objects in H*, dual to the vector fields in TC. As before, the map o 
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depends on the choice of Ti*. Using the L^-pairing with respect to the fixed volume element dx 
and relative to the Euclidean inner product ( • ) in M'^, the momentum map (2.6) is defined for 
images that are functions I & V = by the relation 



Remark 2.1 (Momentum maps). In geometric mechanics, momentum maps generalize the no- 
tions of linear and angular momenta. For a mechanical system, whose configuration space is a 
manifold M acted on by a Lie group G, the momentum map J : T*M g* assigns to each 
element of the phase space T*M a generalized "momentum" in the dual g* of the Lie algebra 
Q of the Lie group G. For example, the momentum map for spatial translations is the linear 
momentum and for rotations it is the angular momentum. 

The importance of the momentum map in geometric mechanics is due to Noether's theorem. 
Noether's theorem states that the generalized momentum J is a constant of motion for the 
system under consideration when its Hamiltonian is invariant under the action of G on T*M. 
This theorem allows us to turn symmetries of the Hamiltonian into conservation laws. 

Remark 2.2 (Momentum of images). Momentum maps for images have been discussed previ- 
ously. In particular, the momentum map for the EPDiff equation of Holm and Marsden [15] 
produces an isomorphism between landmarks (and outlines) for images and singular soliton so- 
lutions of the EPDiff equation. This momentum map was shown in Holm et al. [17] to provide a 
complete parameterization of the landmarks by their canonical positions and momenta. Another 
interpretation of momentum for images in computational anatomy is also discussed in Miller 
et al. [24]. 

We now explain in which sense expression (2.7) is a momentum map. Even though the cost 
functional (2.2) is not invariant under the action of the diffeomorphism group, one may still 
define the momentum map o : V x V* Ti* via 



as done in geometric mechanics, see Marsden and Ratiu [20] and Holm [14]. The action ul is 
defined as ul := dt\t=Q4't ■ / for a curve 0t such that 0o(^) = ^ and dt\t=Q(t>t = u. This is the 
infinitesimal action corresponding to the action of Diff(r2) on V . Although the o-map does not 
provide a conserved quantity of the dynamics, it nevertheless helps our intuition and gives us a 
way to structure the formulas. 

Let us apply this concept to image registration for / G J^(f^), the scalar functions on the 
domain VL. The infinitesimal action is given by 




(2.7) 



(/o/,n) = (/,n/). 



ul = d,l^^{i o 



and thus the momentum map in this case is 



(Jo/,m)v..xV' = (/,-V/-m)= / -{VI-u)fdx={-fVI,u) 



n*xn ■, 



as stated in formula (2.7). The key is to reinterpret the L^-duality between the functions —VI-u 
and / as the duality between the vector fields —fVI and u. 

Using formulas (2.7) and (2.5) in equation (2.4), we regain the stationarity condition (2.3). 
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Remark 2.3. Writing the gradient of the cost functional (2.3) in the geometric form (2.4) has 
several advantages. For example, it allows us to generalize an algorithm that matches images 
as scalar functions, to cope with different data structures, such as densities, vector fields, tensor 
fields and others. Making this generalization allows one to see the underlying common geometrical 
framework in which we may unify the treatment of these various data structures. We can also 
keep the data structure fixed and vary the norm || ■ ||, and thereby alter our criteria of how we 
measure the distance between two objects. 

This geometric framework also enables comparison of different formulations of LDM. For 
example, one may compare the approach from Beg et al. [7] presented here with the symmetric 
approach from Avants et al. [4] and Beg and Khan [6] and the unbiased approach from Hart 
et al. [13], in terms of their respective momentum maps. 

In addition, the geometrical setting introduced here for image analysis allows us not only to 
vary the data structure, but also to change the group of transformations. We will explore this 
possibility in Section 4, when we consider image registration using two diffeomorphism groups 
simultaneously. 

2.2 Abstract Framework 

Diffeomorphic image registration may be formulated abstractly as follows. Consider a vector 
space V of deformable objects on which an inner product ( ■ , ■ ) is defined, that allows us to 
measure distances between two such objects. We can think of V as containing brain MRI images, 
an example frequently encountered in computational anatomy [23]. The distance between two 
objects can be defined as ||/ — = {I — J, I — J) , which in the case of images is the L^-distance 

[ \I{x) - J{x)\'^dx. 
Jq 

The second ingredient is a Lie group G of deformations, that acts on the space V of deformable 
objects from the left 

{g,I) EG xV ^ gl eV. 

In computational anatomy G usually is taken to be the group of diffeomorphisms Diff(f2) or 
variants of it. A diffeomorphism </> G Diff(r2) acts on images by push-forward; that is, by pull 
back by the inverse map, 

(j) . I ■= (f)J = I o or 0-/(x) = J(0-^(x)). 

Roughly speaking, this action corresponds to drawing the image J on a rubber canvas, then 
deforming the canvas by and watching the image being deformed along with the canvas. It is 
also the basis for the familiar Lagrangian representation of fluid dynamics as described in Holm 
et al. [16]. 

Given a curve t ^ gt oi transformations, we define the right-invariant velocity vector G g 

as 

ut = {dtgt)gt^. (2.8) 

We obtain ut by taking the tangent vector of gt and right-translating it back to the tangent space 
at the identity T^G = g, which is the Lie algebra of G. Rewriting (2.8) as 



dt.gt = utgt 



(2.9) 
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and specifying initial conditions at some time t = s, we obtain an ordinary differential equation 
(ODE). If we start with velocity vectors Ut, we can solve this ODE to reconstruct the curve gt. 
This corresponds to the construction of diffeomorphisms as flows of vector fields via the equation 

dt(pt = Uto (j)t, 0o(x) = X. 

Let us denote by g^^ the solution of the ODE (2.9) rewritten as 

dgt,s = utgls, gls = e 

with the initial condition that g^^ is the identity e at time t = s. Since the time t = will play a 
special role, we denote g^ := g^Q. Standard results for differential equations show the following 
properties 

9t,sgs,r = gt.r , gt.s = gtQs^, gt,l = gs,t 

which we will use in our calculations. 

Following the motivation discussed in Section 2.1 we define the abstract version of the cost 
functional (2.2) as 



E{ut) :-- 



i{ut)dt + 



-"g'^lo-hWl 



2ct2 



(2.10) 



where the function £ : g M is a Lagrangian measuring the kinetic energy contained in Ut and 
II ■ II is the norm on V induced by the inner product ( ■ , ■ ). Note that formula (2.10) defines a 
matching problem for any data structure living in a vector space V and any group of deformations 
G acting on V. Although it was inspired by the concrete problem of diffeomorphically matching 
scalar- valued images, the cost function (2.10) no longer contains any reference to image matching. 

Next, we want to deduce (2.4) in our abstract framework. In order to compute the derivative 
DE{ut) we need to know how g"^ behaves under variations 5ut of ut. This is answered by the 
following lemma, the proof of which is adapted from Vialard [33] and Beg et al. [7]. 



Lemma 2.4. Let u : 

Then 



Q, t ^ u{t) be a curve in q and e ^ a variation of this curve. 



d 



Proof. For all e we have 



gls = gls f (Ad,.^ 5u{r)) dr G T,u G. 

=0 Js 



j^gZ = Ueit)gZ. gZ = e. 



Taking the e-derivative of this equality yields the ODE 



d f d 



dt \ de 



gZ 



e=0 



5u{t)gl, + u{t) [ j- 



gZ 



£ = 



and then, using the notation ^^f^^ := ^I^^q^M' compute 



d_ 

dt 



- (gZ) <t)gl {gl) 5gl + fc) (5«(t)<, + u{t)5gl) 

gZ^^(^)9Z 

Adg"^ 6u(t). 
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Now we integrate both sides from s to t and multiply by g^^ from the left to get 

S9ls = 9ls f {Ad,uju{r))dr, 



as required. □ 

Already knowing from (2.4) how the first derivative DE{ut) of the cost functional is going 
to look, we want to establish the necessary notation before we proceed with the rest of the 
calculation. The inner product on V provides a way to identify V with its dual. To / G V" one 
associates the linear form P := {I, ■) &V* . 

Given an action G on V , we define the cotangent lift action of G on vr G via 

(^vr, /) = (tt, g-^l) , for all I eV. 

As mentioned earlier the inverse is necessary to make the dual action into a left action. Finally 
we define the cotangent-lift momentum map o : V ^ V* ^ q via 

(J O TT, m) = (tt, uI) , 

where ul is the infinitesimal action of g on defined by ul = dt\t=ogtI for a curve gt with go = e 
and dt\t=ogt = u. The use of the momentum map was motivated in Remark 2.1. 
Now we are ready to calculate the stationarity condition DE{ut) = 0. 

Theorem 2.5. Given a curve t ut ^ Q, we have 

6£ 

DE{ut) = ^ — (t) = -gPoogliT^. (2.11) 

or, equivalently 

DE{u,) = ^ ^(t) = -1 J° o [gl, (J? - Jl)') , (2.12) 
where the quantities tt, J^, and J} are defined as 

TT := ^ {glh - e V*, jO = g^Io G V, = gl,h G V. 
When G acts by isometrics, the stationarity condition simplifies to 

DE{u,) = ^ ^(t) = -1 J° o (J° - Jiy. 

ou o-^ 

The quantity J^ is the template object moved forward by gt until time t and is the target 
object moved backward in time from 1 to t. 

Proof. Using the notation tt := ^(^'"/o — = ^(^^i ~ •^i)'' ^ may calculate 

{DE{u),6u) = 6 ( f\{u{t))dt + ^ \\g''^Io - hfy 



— {t),6u{t)jdt+U, — 



e=0 
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5i \ 

— {t),6u{t))dt+{TT,6g'^Io) 



Su I 

bu{t)^ dt + (n, (^g- (Ad,.^ ds^ I, 



-(t), 6u{t)^ dt + {{gir' vr, (Ad,.^ 6u{t)) Jo) j dt 

— (t), 5u{t)\ + (jo O TT, Adg.^ 



which must hold for all variations 6u(t). Therefore 

6u 



T-w = -Ad;.^ (/oo(,rrv) 



■gthogl^Ti 



If G acts by isometries, then the action commutes with the flat map and we obtain 



ou a'^ 



□ 



This theorem tells us how to compute the gradient of the cost functional for any data structure 
and any group action. Just like the cost functional (2.10) it is expressed entirely in geometric 
terms and contains no reference to particular examples such as images. This makes the theorem 
widely applicable. 

Remark 2.6. Although the momentum at each time depends on Jq and Ji, it turns out 
that ^{t) obeys a dynamical equation that is independent of /q, /i. The equation in question is 
the Euler-Poincare equation on G. History and applications of the Euler-Poincare equation can 
be found in Holm et al. [16], Marsden and Ratiu [20] and Marsden and Scheurle [21]. 

Lemma 2.7. The momentum ^{t) satisfies 



:W = -<.-W- (2-13) 



- - ad* — ( 
dt 5u 5u 

This is the Euler-Poincare equation on the Lie group G with Lagrangian £ : TG/G ~ g — M. 
Proof. Because the cotangent- lift momentum map is Ad*-invariant we obtain from Theorem 2.5 

^(t) = -g:ioOgl,n 

= -Adl.)-. (/oO«)"V). 
Differentiation of Ad* follows the rules 

dt Ad* ri = Ad* ad* -i ri , 
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9t Ad*-ir/ = - ad* -1 Ad* t] . 

^ 9t ' 9t9t 9t I 

From this we see that 

= ad;Ad; (looig^r'n) 
= -ad„.^(t), 

and so the momentum satisfies the Euler-Poincare equation. □ 

Remark 2.8 (EPDiff equation). When G = Diff(M) the Euler-Poincare equation is the EPDiff 
equation for left action of the diffeomorphisms on the manifold M, 

See Holm and Marsden [15] for a detailed treatment of the EPDiff equation. 

Remark 2.9 (Dependence of lo, h on the initial momentum). It might seem counterintuitive 
that the momentum evolves independently of the objects we are trying to match. However, the 
objects /o,/i do infiuence the momentum |^(t) in a significant way. Namely, solving the Euler- 
Poincare equations requires that we know the initial momentum |^(0) and this initial momentum 
depends on Iq, Ii through the formula 

-(0) = -/oO«)-V. 

Alternatively, we might think of it from the viewpoint of the variational principle. Assume 
that £{u) = -llMp is the squared length of a vector for some inner product ( ■ , ■ ) on g. If we have 



found a vector field Ut and gi, which minimize 

/•I 1 



2 

then the vector field Ut must also minimize 



\V ; 



1 





among all vector fields Ut whose fiows gt coincide with gt at time t = 1, i.e., gi = gi. But this 
means that Ut must be the velocity vector field of a geodesic gt in G. Here we have implicitly 
endowed G with a right-invariant Riemannian metric induced by the inner product ( ■ , ■ ) on g. 
The Euler-Poincare equation (2.14) is just the geodesic equation on the Lie group G with respect 
to this Riemannian metric. 



3 Registration Using the Group of Diffeomorphisms 
3.1 The Setting 

In computational anatomy the group of deformations G is usually the group of diffeomorphisms of 
some domain Q C M.'^. Different types of data used in computational anatomy, such as landmarks. 
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scalar-valued images or vector fields, are deformed by diffeomorphisms via the mathematical 
operations of pull-back and push-forward. Intuitively this corresponds to embedding your data 
into the domain Q, then deforming Q by the diffeomorphism and observing how the data is 
deformed with it. We will go into greater detail about how each of the data types can be 
registered after reviewing some basic notions about the diffeomorphism group. 



Diffeomorphism group 

For technical reasons, we need to consider a group of diffeomorphisms associated to a certain 
Hilbert space of vector fields H. We suppose that His a. subspace of the space of vector fields 
vanishing at the boundary and at infinity, and such that there exists a constant C for which 

\u\i,oo < C\u\n, (3.1) 

where | ■ It^ is the inner product norm of the Hilbert space Ti. and | ■ | is the norm in W^'^{Q). 
Such a Hilbert space defines a unique Kernel K : Q x Q ^ L(M'^,]R'^) such that 



{u,p)l2 = (u, j K{-,y)p{y)dyj 



This also defines a positive, self-adjoint differential operator L (with respect to the L^-inner 
product) such that {u,v)^ = {u,Lv)^2- 

If Ut : [0, 1] — s> 7i is a time-dependent vector field in L^{[0, l],7i), then following Younes [34] 
and Vialard [33], we can consider the solution (pt of the differential equation 

dt(pt{,x) = Uto (ptix), (f)o{x) = X, (3.2) 

and the group 

Gh = {01 I 0i is solution of (3.2) for some Ut E L^{[0, 1],7Y)} . (3.3) 

We shall quickly indicate why Gn is a group, following Trouve [29]. Let 0" and (pi be the flows 
at time t = 1 of the vector fields Ut and Vt- Let Ut '■= —U\-t- Then we have the relation 



0f ° 01 = ri- 



1 1 



since (pl o </)^(x) and (j)i_f-{x) are both integral curves of ut at 0i(x). Taking t = 1, we obtain 
{4>i)~^ = (pi E Gn- To prove that the composition 0" ° 0i is in G?^, we consider the vector field 

/ ^ f 2M2i, if t < 1/2 ^ 

^^^^')-={2d, if t; 1/2' ^^[O'l]- 



In order to compute cp'f'^ , we first solve the ODE for t < 1/2. In this case {u-kv)t = 2u2t ='■ Ut, 
therefore (p^*"" = (p^ = (p2f We then consider the case when t becomes larger than 1/2. In this 

case {u-kv)t = '2v2t~i ='■ Vt and from the situation t < 1/2, we know that at time t = 1/2 the 

-1 



flow 0"*^ takes the value 0". Thus, we must have 0^***' = (p1 o y(p\i2j ° 0i- Now we observe 

that 0j o ^01/2^ ~ 02t-i5 since they are both integral curves of v that coincide at time t = 1/2. 
We thus get the formula 



LU*V IV „ lU 

h =02t-i°0i- 
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Taking t = 1, we get 0^' o = ^^^'^ e Gn- 

Even though Gn is not precisely a Lie group, it comes close enough for our purposes, with TC 
acting as a substitute for the Lie algebra. We can use formal analogies with the finite dimensional 
case to develop applications for computational anatomy. Details about this construction can be 
found in Younes [34], Trouve [29] and results about the regularity of the diffeomorphisms thus 
constructed are found in Trouve and Younes [32] and in Glaunes [11]. 

In the following, when we speak of the group of diffeomorphisms, we will mean the group G-^- 



3.2 Example 1: Landmark Matching 

The simplest kind of objects used in computational anatomy are landmarks. Landmarks are 
labeled collections / = (x^, . . . , x") of points x* G M.'^. Given two sets (x-^, . . . , x"), (y^, . . . , y") 
of landmarks, the landmark matching problem consists of minimizing the energy 



1 1 " 

E{ut) = - \ut\i,dt + —J2m 



i^i-rf- (3.4) 

Our space of deformable objects is ^ = (M"')" with the usual inner product 

n 

(J,J) = 5^x^.y% 

i=l 

for I = (x^, . . . , x"), J = (y^, . . . , y"). The action of the diffeomorphism group G-j-c is by push- 
forward 

0-J:= (0(xi),...,0(x")). 
The corresponding cotangent-lift action on the dual space (M*^")* = M'^" is given by 



and the calculation 



0-/=(D0(xl)-V,...,^0(x'^)-V) 

= ((y\...,y"),(n(xi),...,(x'^ 

n 

= J]/.n(x^) 

i=l 

I n \ 

J^y*5x,,n \ 

«=1 / WxH 



yields the diamond operator 



fx\ 



n 

,x")o(yi,...,y")^ = 5^/5, 

i=l 



where (5x is the delta-distribution defined by / /(y)<5x(y)'^y = /(x) for a test function /(y). 
The condition (2.12) that a minimizing vector field Ut must satisfy is 

1 " 

i=l 
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Consequently, the momentum Lut is concentrated only on the points 0i(x*). By using the Green's 
function K{'k, y) corresponding to the differential operator L, the minimizing condition above 
can be rewritten for the velocity ut as 

n 

= [^0t,i(0i(xO)-^(0i(xO - y')] • 

i=l 

3.3 Example 2: Image Matching 

The large deformation diffeomorphic matching framework used in Beg et al. [7] seeks to match 
two images Iq, Ii by minimizing 

E{ut) = 2 \ut\ndt + ^\\Io ° <PI^ - h\\l2. 

This example has already been discussed in Section 2.1. We review it here by applying the 
abstract formalism developed above. In this example the space V of deformable objects consists 
of real valued functions on Q. We endow this space with the L^-inner product. The group of 
deformations is again the group of diffeomorphisms G^-i, generated by vector fields in Ti. The 
action of Gj^ on V is by push-forward 

for G Gn and I & V. As we have seen, the dual action reads 

■ TT = |detD0"^| (tt o 0^1) 

where |deti3)0| denotes the absolute value of the determinant of D(j). The diamond map in this 
example is 

I on = — ttV/ . 

According to (2.12), a minimizing vector field Ut must satisfy the following necessary condition 

Lut = ^1 det D<l>-l\iJ^ - Jl)WJ^ (3.5) 

where = Iq ° 0j^o, J} = h ° <Pi~i, and (f)t^s is the flow of the vector field Ut 

dt(f)t,s = UtO (f)t^s, 4>s,s{x) = X. 

Equation (3.5) was used in Beg et al. [7] in devising a gradient descent scheme to computationally 
find the minimizing vector field. 

3.4 Example 3: Vector Fields 

Diffusion tensor magnetic resonance imaging measures the anisotropic diffusion of water mole- 
cules in biological tissues, thus enabling us to quantify the structure of the tissue. The measure- 
ment at each voxel is a second order symmetric tensor. It was shown in Pierpaoli et al. [26] and 
Scollan et al. [27] that the alignment of the principal eigenvector of this tensor tends to coincide 
with the fiber orientation in brain and heart. 
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The fiber orientation can be described by a vector field J : f2 — > M*^ and matching two vector 
fields can be formulated as minimizing the energy 

E{ut) = \f \ut\l,dt + ^\\D<P^oI^o<p-^-h\\l,. (3.6) 
Z Jq la 

In this example the space of deformable objects V is the vector space of vector fields in f2, the 
deformation group is the group of diffeomorphisms G-h, generated by vector fields in 7Y, and Gu 
acts on V by push forward 

(J).! = (f)J = L)0o Jo0-i. 

The infinitesimal action of n G on / G V is given by the negative of the Jacobi-Lie bracket 
whose components are 

(97/* ■ f)T^ ■ 

The object dual to vector fields with respect to the L^-pairing are one-forms ir E V* = 
The diamond map is given by 

I OTX = —£i7r — div(J)7r, 



where £i7r denotes the Lie derivative of the one-form vr along the vector field /. In coordinates, 

d 
dx' 

dP ■ d-Ki dP 



writing I = P-^ and vr = iXidx^, we can write the diamond map in the form 



I o-K = — [ 71 j—- + /■'tt— + VTi-— dx\ 
\ ax* ox^ ox^ J 

Using these formulas, we can write the necessary condition for a vector field Ut to minimize (3.6) 
as 

Lut = {£(^t),io + div((0t)*/o)) {\detD(j)tl\{(j)t^i)^7r) , 

where tt = ^ ((0i)*-^o ~ ^ V*- Note that because the b-map does not commute with pull 
backs and push forwards, i.e. 

this formula cannot be significantly simplified. 



3.5 Diffusion Tensor MRI 

Instead of matching only the fiber orientations, we could also match the entire symmetric 2- 
tensor, as was done in Alexander et al. [3] and Cao et al. [9]. In order to do so, we should first 
explain how a diffusion tensor changes under a diffeomorphism. In analogy to images and vector 
fields we could use the push forward by the diffeomorphism. If T is a symmetric tensor-field with 
coordinates Tj^, i.e. 

T{x) = Tijdx' (g) dx^ 

and (p G Gn a diffeomorphism, then the push-forward has the coordinate expression 

0,T(x) = Tij{(f)-\x))Bi{x)Bi{x)dx'' O dx\ (3.7) 

where and Bl{x) is the coordinate matrix of D(f)~^{x). 

In Alexander et al. [3] and Cao et al. [9] a different action was used. At each point a; G ^2 C M"^ 
the orthonormal principal-axis directions ei(a;), 62 (x), 63(0;) of the tensor T(x) are computed, as 



The momentum map representation of images 



December 14, 2009 17 



well as their corresponding eigenvalues Ai(x) > A2(x) > A3(x). Then T can be written as 
T = Aieief + A2e2e^ + Asese^^. The principal axes are each transformed separately as vector 
fields under the diffeomorphisms as in Section 3.4, then normalized and made orthogonal using 
the Gram-Schmidt method. The results are given as: 

^ _ 0*62 - (ei, 0^e2)ei 

110*62 - (ei, 0*e2)ei|| ' 
eg = ei X 

In the above lines, the first principal axis ei is pushed forward by to ei parallel to 0*ei. The 
second principal axis 62 is mapped in such a way, that ei, 62 span the same plane as 0*ei, 0*e2 
and are orthogonal to each other. The last principal axis is then mapped to be orthogonal to 
the first two. The transformed tensor is defined to be: 

■ T = Aieief + A2e2e^ + AsegeJ. (3.8) 

This means, that we transform the principal axis directions as described above, but we do not 
change the eigenvalues. The choice of this action is motivated by the particular application. 
In brain DT-MRI the tensor T{x) describes the diffusivity of water in different directions at a 
position X. The action by diffeomorphisms describes a macroscopic deformation of the brain, 
such as a change of orientation, a growing tumor or a trauma. However, the diffusivity of water is 
governed by the microscopic structure of tissue, which remains unchanged under a macroscopic 
transformation. Therefore, one is looking for a way to transform the tensor, while keeping its 
eigenvalues (the principal diffusivities) unchanged. 

It can be shown that T i-^ 0-T is a left action of Diff (f2) on the vector space of symmetric two- 
tensors. Both of these approaches to Diffusion Tensor MRI given by the actions (3.7) and (3.8) 
have the structure of a Lie group action and thus they may both be cast into our momentum-map 
framework. We leave it to future work to study the different momentum maps that arise for each 
of the actions and the implications that these have for matching of DT-MRIs. 

4 Registration using Semidirect Products 

The examples in the previous section have shown that the abstract formulation of diffeomorphic 
image registration using the diamond operation ( o ) provides a mathematical framework that 
allows us to adapt easily to accommodate different data structures. A second advantage of this 
framework is the ability to perform matching using different groups. The images encountered 
in computational anatomy may contain information on different length scales. Two images can 
vary in their large scale structure as well as in the fine details. In matching such images, it might 
be of advantage to have two groups at our disposal, one to match the large scale behavior and 
the other one to deal with the fine details. This is made possible in our framework by using the 
concept of a semidirect product, which we will review below and then apply in examples. 

4.1 Semidirect Product of Groups 

Consider a Lie group H acting on K from the left by homomorphisms. 

{h,k) e H X K ^ h - k e K, 
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that is, 



hi ■ (/i2 ■ k) = {hih2) ■ k left group action 

h ■ {kik2) = {h ■ ki){h ■ k2) action by group homomorphisms. 



We can then form the semidirect product group G = H (§) K . The group multiphcation in G 
is given by 

9i92 = {h, ki) (/i2, /cs) = {hih2, ki {hi ■ /ca)) 

and the inverse of {h, k) is (/z, k)~^ = h~^ ■ k~^). The Lie algebra g is the semidirect product 
= © ^ of the Lie algebras of H and K. The tangent actions on G are given by 

{hi, ki){h2, k^) = {hih2, ki{hi ■ k2) + ki{hi ■ /cs)) (4.1) 
{hi, ki){h2, k) = [hih2, ki ■ {hik2)) (4.2) 

and the right-trivialization of the tangent bundle is given by 

gg~^ = {h, k){h-\ h'^ ■ k-^) = (jih~\ kk-^ + k{hh-^ ■ k-^' 



The next lemma provides formulas for the adjoint and coadjoint actions of H(§)K on itself and 
its Lie algebra. 

Lemma 4.1. We have the following formulas for the adjoint and coadjoint actions 

Ad^h,k){'^,v) = {AdhV, Adk{h ■ w) + k{AdhV ■ k~^)) (4.3) 

Ad^,,,)(/i, u) = (Ad;(/i + J(fc-V)), . Ad* (4.4) 

eid(vi,wi){v2, W2) = (ad„, V2, ad^„, W2 + vi ■ W2 - V2 ■ wi) (4.5) 

(/^, i^) = (ad*^ fi-wiou, ad^^ u - vi ■ u) , (4.6) 

where J : T*K f)* is the cotangent lift momentum map associated to the action of H on K 

{J{ak),v) = {ak,v- k) , 

and o : X fi* — i> f)* is the cotangent lift momentum map associated to the induced representation 
of H on i 

{w o u, v) := {u, V ■ w) . 

The action {v ,w) E ^ t ^ v ■ w & i is defined as v ■ w = dt\t=Q{h{t) ■ w) for a curve h{t) with 
h{0) = e and dt\t=oh{t) = v. 

Proof. For the adjoint action (Ad) of the group on its Lie algebra, we simply perform the mul- 
tiplications 

Ad(^h,k){v, w) = {h, k){v, w) [h~^, ■ k~^) 
= {hv,k{h-w)) {h~\h-^ -k"^) 
= {hvh-\ k{h ■ w)k-^ + k{hvh-^ ■ k-^)) 
= (Adh V, Adk{h ■ w) + k{Adh v ■ k~^)) 
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and for the coadjoint action (Ad*) on the dual Lie algebra, we pair with (a, 6) G f) (S)6 to define 

(Ad(fc^fc)(/i, iy),{a,b)) = ((/i, z/), Ad(/,,fc)(a, 6)) 

= (/i, Ad,, a) + (z/, Adk{h ■ b) + k{Adh a ■ k'^)) 
= {Adl /i, a) + {h~^ ■ Adl u, b) + (A;~ V, Ad^, a ■ k~^) 
= (Ad* /i, a) + {h'^ ■ Ad* b) + (Ad* (J(A;-V)) , a) 
= { (Ad* (/i + J(A;-V)), /i-i ■ Ad* u) , (a, 6)> . 

For the next identity we differentiate (4.3) and remark that because oi h ■ e = e we get v ■ e = 0. 
Thus, the adjoint action (ad) of the Lie algebra on itself is given by 

d^d(^vi,wi){v2-, W2) = (ad^, V2., ad^^ W2 + Vi ■ W2 + Wi{v2 ■ e) + a.dy^ V2 ■ e - V2 ■ Wi) 
= {ady^ V2, ad^, W2 + Vi ■ W2 - V2 ■ Wi) . 

For the coadjoint action (ad*) of the Lie algebra on its dual, we pair again with {a,b) G ^(§)t to 
see that 

(ad(^i,t«i)(/^' b)) = {{iJ,, u), {ady, a, ad^, b + vi ■ b - a ■ wi)) 

= (ad*j /i, a) + (ad^^ u - vi ■ b) - {wi o u, a) 
= ((ad*^ fi-wiou, ad^^ u - vi ■ u) , (a, b)) 

as stated in the lemma. □ 
If G = H(§)K, the equation 

dtgls = utgls, gls = e- 

can be written as (see (4.1)) 

dt {hi,, kl) = {v,hl, wtkl + vt ■ kl) , hl^ = e, kl, = e, 

where Ut = {vt,Wt) G [)(S)t = g and gf,. = {hl^^k'^.,) G H(§)K. Thus h"^., and kf,. satisfy the 
equations 

dthl = iHhl, dtkl, = wtkl, + ■ kl,. (4.7) 

This means that /i"^ is the flow of the vector fleld vt-, but this is not true for k^, and the vector 
field Wt- The corresponding relation for k^^ is a direct consequence of the noncommutativity of 
the semidirect product. After reviewing these facts about the semidirect product, we will apply 
them to form the semidirect product of two diffeomorphism groups and using this product to 
perform image registration. This is done in the next section. 

4.2 Image Matching with Semidirect Product Groups 

Given a space V of deformable objects, assume that two groups H , K of deformations act on 
V from the left. We imagine H to contain large-scale deformations and K to contain small- 
scale deformations. Since a deformation that captures small structures is also able to capture 
large-scale ones, we will assume that is a subgroup of denoted hy H < K. 

Let us determine the action by group isomorphisms oi H on K subject to the following two 
conditions: 
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• The formula 

{h,k)l:=khl (4.8) 

defines an H (§) K-action on V. Thus h deforms / first on a large scale and then the details 
are captured on a small scale by k. 

• The H(§)K action is effective. If the action is a representation, this means that it is 
faithful. This condition requires that if {h, k)I = I for all I & V, then {h, k) is the identity. 

The first condition implies {hi, ki){h2, k2)I = {hih2,ki{hi ■ k2))I for all hi,h2 G H, ki,k2 G K, 
and I E V. Therefore, kihik2h2l = ki{hi ■ k2)hih2l for all / G V which, by the second 
condition, yields kihik2 = ki{hi ■ k2)hi, that is, the action is necessarily given by conjugation 
[hi ■ /C2) = hik2hi^. In this sense the action by conjugation appears naturally. 

Because of the form of the action on V , the momentum map of the cotangent lifted action of 
H(^K onV xV* has the expression 

/ o TT = (/ Oi TT, J 02 tt) G f)* X r = (f) (§)«)*, (4.9) 

where, I &V , ti &V* , and Joi tt and /02 tt denote the cotangent lift momentum maps of the H 
and K-actions on V , respectively. 

Since the if-momentum map is obtained from the /T-momentum map by restriction, we have 

i*(/ 02 7r) = /oivr, (4.10) 

where I , n ^V* , l : \) ^ t \s the inclusion, and l* -A* ^ P)* is its dual. 
The matching problem using a semidirect product is to minimize the energy 

E{vt,wt)= i{vt,wt)dt+ — WkihJo- hfy, (4.11) 

Jo 2(7^ 

where {hi, ki) are related to {vt, Wt) by 



dtht = vth'l ho = e 

dth = {vt + wt)k'l - k^vt ko = e 



(4.12) 



The last equation is obtained by specializing (4.7) for s = and the action equal to conjugation. 

Theorem 4.2. Given a curve t i— > {vt,Wt) G P)(S)6 the stationarity condition DE{vt,Wt) = for 
the action (4.11) is equivalent to 

-^{t) = -9th Oi 5't,i7r, ^(t) = -gth 02 gt,iir, 
where tt = -^{'gilo — hY and 'gt E K is the solution of the equation 

dtgt = {vt + wt)gt, go = e. (4.13) 

Proof. Let gt = {ht,kt) G H(^K he the solution of the equation dtgt = Utgt, go = e, where 
Ut = {vt,Wt) G [)(S)6 and define gt '■= hht G K. By Theorem 2.5 and (4.9) we get 

— = -g^Io oi gl^-K, — = -g^Io O2 gliix. 
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Since gJo = hhtlo = 9tIo by the definition of tlie if (s) act ion on V, tliis yields 

= -9th oi gt,iT^, ■^{t) = -9th 02 9t,i7^- 
It remains to sliow equation (4.13). By (4.12) we liave 

dtdt = {dtkt)ht + kt{dtht) = {vt + Wt) hht - kiVth + ktVtht = {ik + Wt) gt ■ 
We liave 'qq = k^ho = e. □ 
Tliis tlieorem sliows tliat wlien matcliing witli two groups, tlie momentum j-{t) contain no 



more information than -^{t), since we have 

by (4.10). Nonetheless, this case differs from matching with only one group, since the Euler- 
Poincare equation for the semidirect product reads 

d Si , , Si , ^ Si , , 

which incorporates the Lagrangian of both groups and is genuinely different from the Euler- 
Poincare equation for a single group, which is 

d Si Si 

dts^^'^ = -^<s^^'^- 

4.3 Example: Semidirect Product Image Matching with Two Kernels 

One way of introducing a length scale in image matching is to choose an appropriate kernel 
for the cost of the if-action. If we were to choose for example Lu = u — a^A-u to be the 
differential operator associated to the if^-norm on H, then the corresponding kernel would be 
K(x,y) = e^~'^~^'/"^ where a is a length scale; that is a filter width. A popular alternative 
choice in image registration is the smoother Gaussian kernel K{x,y) = e^"'^"^' \ Increasing 
the value of a increases the cost of forming gradients, or curvature, and thus inhibits nearby 
particles from being deformed differently, while allowing large-scale deformations of the image. 
Sufficiently decreasing the value of a on the other hand would allow fine adjustments in the 
image to be made without requiring much energy cost for the velocity vector field. 

Recall the setting of the example of image matching in Section 3.3. When matching two 
images /q, /i € V := JF(f2) with one kernel, the optimizing vector field Ut satisfies 

u, = ^K* (|detD0,-;i| VJ°(J° - Jl)) , 

where K * f = J K{-, y)f{y)dy denotes convolution with the kernel of the operator L; see (2.3). 

A natural approach for distinguishing between multiple length scales would be to use instead 
the sum of two kernels 



ut = ^ (Ka, + A'aJ * (|det D<P~l\VJ^{J^ - J]] 
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with two length scales ai and a2- We will show, how this approach can be given a geometrical 
interpretation. 

Given two kernel Kai and that correspond to the two length scales ai > a2, we use 
the diagonal Lagrangian i{v,w) = ^\v\'^^ + ll^l^j to measure the energy of the joint velocity- 
vector {v,w). The norm | ■ |q,- is associated to the inner product coming from the kernel Ka^, 
i = 1,2. We assume that the associated Hilbert spaces Hai C Ti^^ verify the hypothesis (3.1). 
Let C be the groups associated to Tia^, via (3.3). The element {ip,ri) G Ga^ (§)Ga2 
acts on V = JF(f2) by the action (4.8), that is, 

{ip, ri) ■ I := [j] o if)) ■ I = I o (rj o ip)~^ = I o ip~^ o r]~^. 

The matching problem with the semidirect product group Gcei (§) Ga2 is to minimize the energy 

E{vt,wt) = - \vt\l, + \wt\l^dt + — \\lo o ip-^ o 7]^^ - h\ 



|2 



By Theorem 4.2, the energy is minimal if 

Vt = K^^ * (-4>th Oi 0t,i7r) , Wt = Ka2 * (^-4>th 02 0t,ivr) 

and ^ ^ 

5t0t = {vt + Wt) o (j)t, 00 = id. 

The example of single kernel image matching in Section 3.3 showed us that 

-4>tho4>t,n^ = ^|detD0,7i^|VJ?(J? - Jl), 

with J° = Jo o 0~o\ J} = 

^1 ° 0t 1 • -By denoting ut :— Vt + Wt the velocity vector field of (pt, we 

see that 

ut = ^ {K^, + K^2) * (idet D4>il\S/JKJt - Jl)) ■ (4-14) 

Thus, matching images with the sum of two kernels corresponds to using a semidirect product 
of diffeomorphism groups. This provides a geometrical interpretation for an approach that was 
suggested intuitively. 



5 Symmetric Formulations of Image Registration 

The cost functional (2.10) is not the only choice possible in the large diffeomorphism matching 
framework. Other cost functionals have been proposed in the literature, which make the regis- 
tration problem symmetric. A consequence of the choice (2.10) is that it matters, whether we 
choose to register Iq to Ji or vice versa. In some applications it may be useful to conceptually dis- 
tinguish between Jg and Ji, for example when the template Jq is available in a higher resolution, 
but in other cases one would prefer a symmetric cost functional to (2.10). Such cost functionals 
have been proposed in Beg and Khan [6], Avants et al. [4] and Hart et al. [13]. We will show 
how they can be analyzed geometrically similar to the cost functional (2.10) in Section 2. 

Example 5.1. The approach described in Avants et al. [4] and Beg and Khan [6] can be ab- 
stractly described with the following cost functional 

/■I 1 

E{ut)= l{ut)dt+—\\gao-g^_^Ji\\l 
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where gt^s is the flow of ut. Since we now evaluate the inexactness of the matching in the midpoint 
i = i of the interval, this choice of the cost functional leads to a symmetric formulation of LDM. 
A calculation similar to theorem 2.5 with vr = -^{gilo — gi ihY 



{DE{ut) 



5ut) = (^(^)' ^"(*)) + - Sgiji) 



(^W'^^W)^^ 



71, gi i^j^ Ad^-i 5u{t)dt^ Jo - ^1^1 i^j^ Mg,, 5u{t)dt^ h 



— {t),6u{t))dt 
du I 



+ / (lQog^'^Ti,Mg-i5u{t)^dt+ j (lio g^^^'K,Mg^ ,5u{t)^ dt 
Jo J ^ 

= (^^{t)+gtIoOg^iTi,6u{t)'^dt + (^^{t) + gt^ihog^iTx,5u{t)^dt 

shows that a minimizing vector field has to satisfy 

^{t) = -gthog,.Tx, te [0,1/2] 
— {t) = -gt,ihog^ i7r, t G [1/2, 1] 
n = j-,{giJo-9i^,ihf 

The momentum has a structure very similar to that of Theorem 2.10, but now we have a discon- 
tinuity at time t = 1/2. 

Example 5.2. Another approach to symmetrize the registration problem considered in Beg and 
Khan [6] is via the cost functional 

1 1 .1 



E{ut)= J £{ut)dt+— J Wgth - gt,ih\\vdt 

Instead of minimizing the matching error at some chosen time like t = for the classical LDM or 
t = I as in the previous example, this approach averages the error over the whole time-interval. 
Using the notation vrt = -^{gtlo — 9t,ilif we can again calculate the derivative of E{ut) 



{DE{ut),6ut) 



Jo (^*^^'*'^^*^^'') ^ / {nr,SgrIo- Sgr,ili)dr 



(^W'^^W)^^ 



+ (nr,gr (^J^ Mg-1 5u{t)dt^ Iq - gr^i (^J^ Ad^,, 6u{t)dt^ h 



(^W'^^W)^^ 
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pi pr pi pi 

+ / {gJoO 9t,rT^r,^u{t))dtdr + / / {gt^JiO gt^rT^r,^u{t))dtdr 

Jo Jo Jo Jr 

+ / {gJoO 9t,r'n-r,Su{t))drdt + / / {gt^iho gt^rTTr,Su{t))drdt 
Jo Ji Jo Jo 

^(^) + J {9t,ih'i-[o,t]{r) + gJoMtA]i^))'^9t,r'^rdr,5u{t))dt. 
This calculation yields necessary conditions for the minimizing vector field 

j^{t) = - J {9t,ih'i-[o,t]ir) + gtIolit,i]ir)) o gt,rTirdr , 



1 



= -^{.9tIo - 9t,ih) 



Here, l[o,j](r) is the indicator function of the interval [0,t], i.e. l[o,t](r) = 1 for r e [0,t] and 
otherwise. The momentum again displays the momentum map structure, in this case involves 
an average over time. 

Example 5.3. A third approach to symmetric registration was presented in Hart et al. [13]. 
They proposed to allow inexactness in the initial as well as in the final image by choosing the 
cost functional ^ 

E{uu I) = i{ut)dt loWl + ^\\9il - /illy- 

This cost functional treats J G as an additional free variable. Intuitively this approach means 
that we are looking for an energy minimal path such that both the starting and the ending points 
match Jq and Ji as well as possible. Computing the necessary conditions for the pair {ut, I) to 
be minimizing with vr = -^{gil — yields 

/ 5i \ 1 

{DE{uuI),i6ut,5I)) = J (^ — it),6u{t))dt+^{l'-ll5I) + {n,5gJ + gi5I) 

= dt + ^(/' - /o + ^'9i'^, SI) + 91 {^f^ Ad^-i 5u{t)di)j I 

= (^ — {t)+gtIogt,in,5u{t)\dt+—{I^-I^, + a^g^^7r,5I). 
This leads to 

-^{t) = -9tIo9unr, 
n = ^igj-hf. 

If we are dealing with images I G JF(f2,]R) as in section 3.3, then the equation for can be 
solved explicitly to find 

_ /q + |-D0i|/i o 01 

i + m^ ■ 

In this case / constitutes a weighted average of Iq and the deformed image ■ Ji at time t = 0. 
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We have seen in all these examples that the momentum had a very similar structure 
based on the momentum map. They differed in the time point at which the inexactness of the 
matching was measured, or, as in the last case, in which image was being compared. We have 
restricted our attention primarily to only one of these possible formulations of LDM. However, 
the geometric interpretations are similar in all cases and the momentum map plays the dominant 
role in each case. 



6 Nonlinear Generalizations 

We now show that the formalism developed in §2.2 generalizes easily to the case when the set 
of images is not necessarily a vector space and the cost function is not necessarily the Euclidean 
distance. This situation arises, for example, in the Landmark Matching Problem associated to 
points on the sphere for the study of neocortex, see Miller et al. [23] and references therein. 

Suppose the set of images is a manifold Q on which a group of transformation G acts on the 
left. As before, we denote by gl the action g E G on I E Q. We consider a cost function of the 
form ^ 

E{ut)= ! t{ut)dt + F{glh,h), (6.1) 



where F is defined on Q x Q. When Q is a vector space V with inner product norm || ■ ||y, we 
recover the cost function (2.10) by choosing 

The next theorem establishes the stationarity condition associated to the cost in (6.1). 
Theorem 6.1. Given a curve t ^ Ut in the Lie algebra g of G, we have 

DE{m) = ^ -{t) = - J {gl, d,F{Jl h)) , 

where J : T*Q g* is the cotangent bundle momentum map and diF(J?,Ii) G T*oQ is the 

■■'i 

tangent map to F relative to the first variable. The momentum ^{t) satisfies the Euler-Poincare 
equation 

dt6^^'^ = -'^-^6^^'^- 

Proof. The proof is similar to that of Theorem 2.5. We will use the formula {3{aq),u) = 
{ag,UQ{q)) for the momentum map J : T*Q g* associated to the cotangent lift action. Using 
Lemma 2.4, we calculate 



{DE{ut), 5ut) = 5 (^^' l{u{t))dt + F ((^r/o, h\ 



^(t), 5u{t) ) dt + {d,F{Jl h), i6g'i)Io) 



Jo [{^(t)^^''(t))dt+(^i9ir'd,F{Jlh), (Ad,./n(t))^/o^^ dt 
I' + m-'^^F{J'„h)) , Ad,. Ju{t))^ dt 
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which must hold for all variations 6u{t). Therefore, 

|l(t) = -Ad;„^ (J {{g-r'd,F{j'„h))) 
= -3{gl,d,F{Jlh)) , 

as required. The same proof as for Lemma 2.7 shows that the Euler-Poincare equations are 
verified. □ 

When Q is a vector space V , this stationarity condition can be rewritten by using the o map, 

as 

^^{t) = -j',o{gl,d,F{Jlh)). 

Landmark matching on manifolds. In the case of the Landmark Matching Problem on a 
Riemannian manifold Q, one chooses the cost function 

n ^ 

F{qi, ..,qn]Pi, ...,Pn) ■■= ^ 7r^d{qi, p,^'^ 



2a2 

i=l 



I J 1 



where d is the Riemannian distance. This approach is used for imaging of the neocortex, where 
Q is taken to be the sphere S"^. The energy to minimize has the form 



1 " 1 

E{ut) = - \ut\n + ^^d{M(li),Pi) 

where qi,Pi C S"^ are given. 



1=1 



LDM multimodal image matching. The framework developed above allows us to under- 
stand geometrically the model developed in Vialard [33], §3.2. This model deals also with a 
change of intensity in the image I : Q ^ X. This change of intensity can be modeled by an 
action 77 o / of a diffeomorphism of the template co-domain X. In this case, the energy can have 
the general form 

E{vt,wt) = / £{vt,wt)dt + F{r]io Iqo (f)~^,I^), 
Jo 

where rjt G Diff(X) and (pt G Diff(f2) are the flows of Vt and Wt, respectively. This problem can 
be recast in our formulation by considering the action of the direct product Diff(f2) x Diff(X) 
on the manifold Q = J-'{Q, X) given by 

(0, 77) ■ / := ?7 o / o 

For simplicity, we suppose that X is a vector space, but in general X can be an arbitrary manifold. 
The cotangent lifted action on tt reads 

(0, f]) ■ (/, 7r)(x) = I det Dr\x)\Dr]~\l{4>{x)) f ■ 7r{4>~\x)) 
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and the momentum map is 

J(/, tt) = ^-vr • V/, j 7r(x)5/(^)dx^ . 
Using these formulas, the stationarity condition is 

where J° = r]t o Iq o (f)~^ and 

Mx) ■.= \detD<P~lix)\ {Dr^tAJK^))y''d,F{j',,I,){<p-l{x)). 
The last expression is obtained using the formula of the cotangent lifted action and the equality 

For more discussion, see Vialard [33]. 

Alternative approach. We now consider an alternative approach that affects the geometric 
shape of the image / : — > X, as considered in Trouve [29]. This approach is different from that 
considered above. For example we can consider the case X = S*^ of images of unitary vectors 
m R^. In this case the shape can be modified by letting various groups of matrices acting on 
5*^. These matrices are of course allowed to depend on the domain Q. We thus need to consider 
the group J^{Q,G), where G is a group acting on X. In order to also take into account the 
transformation on the domain, the semidirect product DiS{Q)(§)J-'{Q,G) 3 (0,6') needs to be 
considered as in Trouve [29]. This group acts in a natural way on the space J-'{Q,X) of images 
via the left action 

(0,^)-/= (e/)o0-l, 

where the function 61 is defined by {61){x) := 9{x)I{x) and in the last term we use the G-action 
on X. A vector field on this Lie algebra has components (m, v) where m is a vector field on Vt and 
p : VL ^ Q. Using the multiplication rule (0, 6')(0, 9) = {(po 0, (6^ o 0)^) in the semidirect product, 
the ODE dt{(j)uOt) = {ut,iyt){(j)t,Ot) reads 

0i = Ui o 04, 9t = {vt o (t)t)Ou 00 = e, 9o = e. 

For simplicity, we suppose that A is a vector space. The infinitesimal action on the space of 
images reads (u, v)I = vl — V/ • u hence the cotangent bundle momentum map is 

J(J,7r) = (-TT • V/,/o7r) , 

where /ovr is the function with values in g* defined by (/o7r)(x) = I{x)oti{x) and the diamond on 
the right denotes the momentum map associated to the action of G on A. In order to formulate 
the stationarity condition, we also need the expression of the cotangent lifted action given by 

(0, e) ■ (J, vr) = {{91) o 0-1, 1 det D(t)~^\ {Otx) o 0-i) . 

The cost function has the form 



E{ut,yt)= [ i{ut,ut)dt-F{{9Jo)o(p-\l,). 
Jo 
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The stationarity conditions are thus given by 



OU Of 



where = 6t) ■ Jq = (OJo) ° (pt ^ and 



detD<P]:i\ {9i,tdiF{Jlh)) 



For example, when F{I, J) = — J\\l'^, relative to an inner product on X, then diF{I, J) = 

i(/ - Jf e J^{n,X*), where b is associated to the inner product on X. In this case, the 
stationarity conditions are 

^{t) = I det D<j^^}\ ( J° - J^fVJ?, ^ = J?o I det D<j^^}\ ( J° - J^f. 

OU ' Of ' 



7 Conclusions 

This paper has revealed that Beg's algorithm from Beg [5] and Beg et al. [7] for image regis- 
tration in the LDM framework is the cotangent-lift momentum map associated to the action of 
diffeomorphisms on scalar functions. Accordingly, the momentum map has emerged as a central 
organizing principle in the abstract framework inspired by image registration. The momentum 
map provides the means of unifying the LDM approach for the registration of different data 
structures that use different penalty terms and different Lie groups. Different data structures 
summon different group actions to define their transformations and they will therefore give rise 
to different momentum maps. But once the momentum map is computed, it is straight-forward 
to implement the corresponding gradient-descent scheme for image registration. Although not 
emphasized in this paper, it should be noted that the specification of distance on the space of 
images will also be reflected in the momentum map. Exploring the specification of distance and 
dealing with other data structures has been left for future work. For example, the pioneering 
work of Alexander et al. [3] and Cao et al. [9] on the registration of DT-MRIs led to the action 
on symmetric tensors discussed in Section 3.5. We plan to compare the momentum map for 
this action with the usual push-forward action on tensor fields to gain further insights into the 
matching procedures for tensor data structures. 

Images encountered in applications often contain information at several length scales. A 
heuristic approach for adapting the registration procedure to take into account these length scales 
suggested replacing the kernel in (2.3) by the sum of two kernels K^^ + K^^i with two different 
length scales ai and 02 for their corresponding filters. We have shown that this strategy has 
a geometric interpretation. Namely, instead of using a single diffeomorphism group to perform 
image registration, we can use the semidirect product of two such groups, each associated to its 
own length scale. The resulting equations (4.14) then coincide with the sum-of-kernels strategy. 
Similarly, the same result could be obtained for the sum of three and more kernels. 

Other formulations of LDM that were intended to make the registration symmetric, as pro- 
posed by Avants et al. [4], Beg and Khan [6] and Hart et al. [13], were also written geometrically. 
We have shown that all these cases exhibit similar momentum map structures. The main differ- 
ence arises from the choice of the time at which the momentum map is to be evaluated. Once 
again, the momentum map appears as a unifying framework allowing systematic comparisons 
among the different examples. 
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We have also explored a natural generalization of the framework to incorporate data structures 
living in manifolds, which do not have the linear structures of vector fields. Examples included 
landmarks on a sphere. Since in this case no norm is available to measure distances between 
two images, a distance function must be chosen. Further applications and capabilities of this 
nonlinear framework will be explored in future work. 
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